$title  CF_CO2_NH.GMS

* ###################################

* Per capita income, consumption patterns, and CO2 emissions

* Journal of the Association of Environmental and Resource Economists

* Justin Caron and Thibault Fally

* April 2021

*###################################

* SIMULATION FILE with NH CES preferences (CLM)

*###################################


$set SLASH \
$if %system.filesys% == UNIX $set SLASH /


* declare set of demand estimates to use : theta4, notc, tc

$if not set demandest $set demandest theta4
* need to choose H or NH here

* options: NH_CLMstd  
* options: NH_CLMquad 
* options: NH_CLMshift

$if not set spec $set spec NH_CLMshift

* declare supply side specification: noresprod, resprod, resprod_06, resprod_15,  resprod_50
*$if not set prodspec $set  prodspec resprod_06
*$if not set prodspec $set  prodspec resprod_15
$if not set prodspec $set  prodspec resprod_50


$if not set ds $set ds gtap8gas



* TWO COUNTERFACTUALS: CHOOSE ONE:
parameters income_shock, tcost_shock, income_shock_toresourcefactor;
tcost_shock = 1;
income_shock = 1.01;

* set this to 1 if productivity of resource factor is to be shocked as well -- 0 if off
income_shock_toresourcefactor = 1;
$if "%prodspec%"=="noresprod"  income_shock_toresourcefactor = 0;

* CHOOSE NB OF ITERATIONS:  (use 20 for final results, though converges pretty well after ~5)
* note: real variables (in quantities) converge. Nominal variables not necessarily.
Sets
T iteration  /T1*T20/;

Sets
r country  /
AUS,NZL,CHN,HKG,JPN,KOR,MNG,TWN,KHM,IDN,LAO,MYS,PHL,SGP,THA,VNM,BGD,IND,NPL,PAK,LKA,CAN,
USA,MEX,ARG,BOL,BRA,CHL,COL,ECU,PRY,PER,URY,VEN,CRI,GTM,HND,NIC,PAN,SLV,AUT,BEL,CYP,CZE,
DNK,EST,FIN,FRA,DEU,GRC,HUN,IRL,ITA,LVA,LTU,LUX,MLT,NLD,POL,PRT,SVK,SVN,ESP,SWE,GBR,CHE,
NOR,ALB,BGR,BLR,HRV,ROU,RUS,UKR,KAZ,KGZ,ARM,AZE,GEO,BHR,IRN,KWT,QAT,SAU,TUR,ARE,EGY,MAR,
TUN,CMR,CIV,GHA,NGA,SEN,ETH,KEN,MDG,MWI,MUS,MOZ,TZA,UGA,ZMB,ZWE,BWA,NAM,ZAF,ISR,OMN
/;

alias (r,s);
alias (r,rloop);

Sets
i industry  /
isr,obs,ros,osg,pdr,wht,gro,v_f,osd,c_b,pfb,ocr,ctl,oap,rmk,wol,frs,fsh,coa,oil,gas,omn,cmt,omt,vol,mil,pcr,
sgr,ofd,b_t,tex,wap,lea,lum,ppp,p_c,crp,nmm,i_s,nfm,fmp,mvh,otn,ele,ome,omf,ely,wtr,cns,trd,otp,wtp,atp,cmn,ofi
/;

alias (i,j);

Sets
f factor /Land,UnskLab,SkLab,Capital,NatlRes/;


* Loading data:
parameter fittedexp, coeffs, U, scale, sigma, theta, beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop, incelast;

$if "%spec%"=="NH_CLMshift"  $gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_CLM_flex_shifter.gdx
$if "%spec%"=="NH_CLMquad"  $gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_CLM_flex_quadratic.gdx
$if "%spec%"=="NH_CLM"  $gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_CLM.gdx

$load  coeffs  fittedexp

$if "%spec%"=="NH_CLMshift"  $load U
$if "%spec%"=="NH_CLMquad"   $load U
$if "%spec%"=="NH_CLM"   U(r) =1; 


* And everything else from DATAPREPARATION.GMS
$gdxin estimates%SLASH%co2_data_tc_avg.gdx
$load  beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop, incelast

theta(i) = 4;

* -------
* Theta 4:  H, sigma = scale, NH, sigma = sigma * scale
*  here, the average sigma level is recoverd from Mu

* TC: H, we don't know the average level of sigma
* back it out assuming the average level of theta = 4

$if "%demandest%"=="tc" theta(i)$coeffs("backed out theta","coeff",i) = coeffs("backed out theta","coeff",i);

*normalize theta's such that the average equals = 4 and min = 1:
theta(i)$(theta(i) = 0)  = 0;
theta(i)$(theta(i) < 0)  = 0;
theta(i) = 4 * theta(i) * sum(j,1) / sum(j,theta(j));
theta(i)$(theta(i) = 1)  = 1;
theta(i)$(theta(i) < 1)  = 1;
theta(i) = 4 * theta(i) * sum(j,1) / sum(j,theta(j));

* some theta's are very large, set max =12
theta(i)$(theta(i) > 12)  = 12;
* only affects 1 sector


* load clm parameters
parameter rho, b, logUsquare;

* for sigma pick one, they are all the same
sigma =  coeffs("sigma","coeff","obs");

* note: estimated rho is actually rho *(1-sigma)  
* in model we need just rho so need to transform
$if "%spec%"=="NH_CLMquad" rho(i) =  coeffs("rho","coeff",i) / (1-sigma);

* note: b = nu and was estimated with minus sign
$if "%spec%"=="NH_CLMquad" b(i) =  coeffs("nu","coeff",i) / (1-sigma);

*** NOT THE CASE IN THE SHIFTER SPEC
$if "%spec%"=="NH_CLMshift"  rho(i) =  coeffs("rho","coeff",i);
$if "%spec%"=="NH_CLMshift"  b(i) =  coeffs("b","coeff",i);


* standard CLM:

* note: estimated rho is actually (eps-sigma)
* in model we need eps-sigma

* note: shoule be negative as it corresponds to (eps-sigm)/(1-sigma)
$if "%spec%"=="NH_CLM" rho(i) =  coeffs("rho","coeff",i) ;

* U is loaded

* for quadratic form, pre-compute logUsquare
logUsquare(r) = (log( U(r)))*(log( U(r)));

* check that rho is positive (mostly)
display U, logUsquare, rho;


set not_fd(i)  sectors with no sigma estimate (where dropped from estimation);
not_fd(i)$(not rho(i))  = yes;

* set theta = 4 for sectors for which we have no estimates
theta(i)$(not theta(i)) = 4;

display theta, sigma;


* ----
* CHOICE OF PRODUCTION SHARES 

* compute fitted production shares from fitted demand shares
* idea: keep observed values for all non-demand related values


parameter finaldemand;

finaldemand("true",i,r) = fd(i,r); 
finaldemand("non-homoth",i,r) = fittedexp("non-homoth",i,r); 
finaldemand("homoth",i,r) = fittedexp("homoth",i,r); 
 

** REDEFINING ABSORPTION AND PRODUCTION FOR CONSISTENCY:
parameters production, absorption;

production(i,r) = sum(s, trade("data",i,r,s));
* INSTEAD OF:	 
*production(i,r) = vom(i,r);
		 
absorption(j,r) = fd(j,r) + sum(i,intersh(j,i,r) * production(i,r));
* INSTEAD OF:
*absorption(i,r) = sum(j, vdfm(i,j,r)+ vifm(i,j,r));
* THIS GUARANTEES WELL DEFINED GOSH COEFFICIENTS -- BUT CONSISTENCY ISSUES WITH IMPORT FLOWS.


display absorption, production;
* Production: some zeros
* Absorption: zeros with previous definition


parameter
dem_sh(i,r)  share of sector i in total final demand
trade_sh(i,r,s)  share of country r in country s's absorption (Import shares)
prod_sh(i,r,s)  share of country s in country r's sales
fact_sh(f,i,r)   share of sector i in demand for factor f
endow_sh(f,r) share of factor f in total endowment (income)
gosh_m(j,i,r) gosh matrix coefficient (how much j is used in industry i)
gosh_f(j,r) gosh matrix coefficient (how much j is used by final consumers);


parameter prod_sh_adj;

parameter vom_adj , vfm_adj;


* Loading FITTED PRODUCTION SHARES: from CF_compute_fittedprod_shares.gms

set fit /non-homoth, homoth, true/;
parameter fitted_prod_sh(i,r,s,fit), fitted_production_absorption;
$gdxin estimates%SLASH%fitted_prod_shares_%spec%_%demandest%.gdx
$load  fitted_prod_sh fitted_production_absorption

* note: HOMOTHETIC not done but would be easy

* here, also wand to adjust demand shares, employement shares, and gosh coeffs

** SPECIFICATION:**** chose here:
* 1 - bilateral production shares
*prod_sh_adj(i,r,s) = fitted_prod_sh(i,r,s,"true");
prod_sh_adj(i,r,s) = fitted_prod_sh(i,r,s,"non-homoth");
*  tested that using only fitted production shares works: eq solutions is 1

* 2 - demand shares
*dem_sh(i,r) =  finaldemand("true",i,r)  / sum(i.local, finaldemand("true",i,r) );
dem_sh(i,r) =  finaldemand("non-homoth",i,r)  / sum(i.local, finaldemand("non-homoth",i,r) );


* 3- absorption and production  --> WILL BE USED TO COMPUTE GOSH COEFFICIENTS 
* keeping observed:
*production(i,r) = production(i,r);
*absorption(j,r) = absorption(j,r);

* or better: fitted (prob better but didn't try)
*production(i,r) = fitted_production_absorption("production","true",i,r);
*absorption(j,r) = fitted_production_absorption("absorption","true",j,r);

* or NH:
production(i,r) = fitted_production_absorption("production","non-homoth",i,r);
absorption(j,r) = fitted_production_absorption("absorption","non-homoth",j,r);

* 4- FACTOR SHARES: also require production and absorption to be computed
* keeping observed (use this to minimize distortions)
*vom_adj(i,r) = vom(i,r);
*vfm_adj(f,i,r) = vfm(f,i,r);

* using 'fitted' absorption/production (from true demand)  (prob better but didn't try)
*vom_adj(i,r) = fitted_production_absorption("production","true",i,r);
*vfm_adj(f,i,r) =  beta(f,i, r) * vom_adj(i,r);


* or NH
vom_adj(i,r) =  fitted_production_absorption("production","non-homoth",i,r);
vfm_adj(f,i,r) = beta(f,i, r) * vom_adj(i,r);

* -------------------------------------------------------------------------

parameters CHECK_prod_sh;
CHECK_prod_sh(i,r,fit) = round(sum(s,fitted_prod_sh(i,r,s,fit))-1,13);

display CHECK_prod_sh;


* j being used in industry i:
gosh_m(j,i,r)$absorption(j,r) = intersh(j,i,r) * production(i,r) / absorption(j,r);

gosh_f(j,r) = 1 - sum(i,gosh_m(j,i,r));
* normalized to unity when absorption is zero


display intersh, gosh_m, gosh_f;
* Some negative values because of rounding errors, not a problem 

* defined as: from r to s:
trade_sh(i,r,s)$sum(r.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(r.local , trade("data",i,r,s));

parameter comp_ab1, comp_ab2;
comp_ab1(i,s)$absorption(i,s) = round(sum(r,trade_sh(i,r,s)) - 1,14);
comp_ab2(i,s)$(round(sum(r,trade_sh(i,r,s)) - 1,14)) = absorption(i,s);
display comp_ab1, comp_ab2;

trade_sh(i,s,s)$(sum(r.local,trade("data",i,r,s))=0)  = 1;

* this is to be replaced if using fitted shares
* defined as: from r to s:
prod_sh(i,r,s)$sum(s.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(s.local , trade("data",i,r,s));

* adjustment when sum is zero:
prod_sh(i,r,r)$(sum(s.local , trade("data",i,r,s))=0)  = 1;

fact_sh(f,i,r) =  beta(f,i, r) * vom_adj(i,r) / sum(i.local,  beta(f,i,r) * vom_adj(i,r))  ;

endow_sh(f,r) = sum(i,vfm_adj(f,i,r)) / sum((f.local,i),vfm_adj(f,i,r)) ;

* CHECKS TO ENSURE THAT CHG=1 IS A SOLUTION:

parameters CHECK_budget, CHECK_absorption, CHECK_production, CHECK_supplier, CHECK_index, CHECK_fmc, CHECK_income;
CHECK_budget(r)= round(sum(i,dem_sh(i,r)) - 1,14);
CHECK_absorption(j,r) = round(gosh_f(j,r) + sum(i$gosh_m(j,i,r),gosh_m(j,i,r)) - 1,14);
CHECK_production(i,r) = round(sum(s,prod_sh(i,r,s)) - 1,14);
CHECK_supplier(i,r) = round(prod(f$beta(f,i,r),1**beta(f,i,r))**theta(i) * prod(j,1**intersh(j,i,r))**theta(i) - 1,14);
CHECK_index(i,s) = round(sum(r,trade_sh(i,r,s)) - 1,14);
CHECK_fmc(f,r) = round(sum(i,fact_sh(f,i,r)) - 1,14);
CHECK_income(r) = round(sum(f,endow_sh(f,r)) - 1,14);

display CHECK_budget, CHECK_absorption, CHECK_production, CHECK_supplier, CHECK_index, CHECK_fmc, CHECK_income;
* ISSUES FOR CHECK_production AND CHECK_index
* All fine 

parameter trade_sh_adj;
trade_sh_adj(i,r,s)$absorption(i,s) = trade("data",i,r,s)/absorption(i,s);


***---------------------------------------------------------------------------------------------------

* ADJUSTMENT TO GET A CES PROD FUNCTION IN NATURAL RESOURCES:
* may 2017

Set	f_noRes(f)	 factor other THAN NatlRes
	res_sect(i)   sector with sector specfic resource 
	pe(i)		 primary energy 
		;

f_noRes(f) = yes;

f_noRes("NatlRes")  = no;

pe("gas") = yes;
pe("oil") = yes;
pe("coa") = yes;


* ADJUSTMENTS 1 -- for sectors to take out the natural resource factor for non-enegy sectors INCLUDES OMN FSH FRS
* transfer that resource to capital
* i.e, only keep NATLRES for the ENE sectors

* solutions:
* FLAG
* 1) make a new factor, tradable across these sectors..  pro: simple   con: not realistic
* 2) merge with other factor.. FOR NOW THIS -- capital
* 3) keep sector specific, but with supply elasticity = 1 ..  Problematic, we don't really want to fool with those sector's behavior


set nonatlres(i) / OMN, FSH, FRS /;


* change vfm and beta
vfm_adj("Capital",nonatlres,r) = vfm_adj("Capital",nonatlres,r) + vfm_adj("NatlRes",nonatlres,r) ;
beta("Capital",nonatlres,r) = beta("Capital",nonatlres,r) + beta("NatlRes",nonatlres,r) ;
vfm_adj("NatlRes",nonatlres,r) = 0;
 beta("NatlRes",nonatlres,r) = 0;
 
* re-compute factor and endowment shares
fact_sh(f,i,r)$sum(i.local,  beta(f,i,r) * vom_adj(i,r))  =  beta(f,i, r) * vom_adj(i,r) / sum(i.local,  beta(f,i,r) * vom_adj(i,r))  ;
endow_sh(f,r) = sum(i,vfm_adj(f,i,r)) / sum((f.local,i),vfm_adj(f,i,r)) ;


* ADJUSTMENTS 2 -- some countries do not produce certain fossil fuels at all.
* in this case, best to fix W_res to 1
set no_ene_prod(i,r);

no_ene_prod(pe,r)$(not production(pe,r)) = yes;


parameter beta_all(i,r), beta_Res(i,r), beta_fnores(f,i,r), endow_sh_res  ;

* note: as of now, the beta are NOT COUNTRY-SPECIFIC

beta_all(i,r) = sum(f, beta(f,i,r)) ;

beta_fnores(f,i,r) = beta(f,i,r) / sum(f.local$f_noRes(f), beta(f,i,r));

beta_Res(i,r) = beta("NatlRes",i,r) / sum(f.local, beta(f,i,r));




parameter	supply_elast(i,r)    targeted elasticity of supply
		nu(i,r)  elasticity of subsitution between resource and non-resource factor  -- to be calibrated to reach desired  supply elasticity ;


supply_elast("oil",r) = 0.75 ;

* necessary to keep nu positive - cannot go much below this
supply_elast("gas",r) = 0.75;
supply_elast("coa",r) = 0.8;


* minimum value for coal to have nu>0 is 0.8. set other two to 0.75
$if "%prodspec%"=="resprod_06" supply_elast(pe,r) = 0.75 ; supply_elast("coa",r) = 0.8;
$if "%prodspec%"=="resprod_15" supply_elast(pe,r) = 1.5;
$if "%prodspec%"=="resprod_50" supply_elast(pe,r) = 5;


* formula from the appendix
nu(pe,r) = (supply_elast(pe,r) +1 - 1/beta_all(pe,r)) * (beta_all(pe,r) *  beta_Res(pe,r)    /  (1 - beta_Res(pe,r))) ;

res_sect(i)$(sum(r, beta_Res(i,r)))  = yes;

* for calculation of change in income

endow_sh_res(i,r)$res_sect(i) = vfm_adj("NatlRes",i,r) / sum((f.local,i.local),vfm_adj(f,i,r)) ;

* check income again

CHECK_income(r) = (sum(f$f_noRes(f),endow_sh(f,r)) + sum(i$res_sect(i),  endow_sh_res(i,r) ));

***---------------------------------------------------------------------------------------------------

parameter bitc_shock;
* Trade cost shock just between different countries:
*bitc_shock(r,s) = tcost_shock;
bitc_shock(r,s) = 1;


positive variables
                chg_income(r) per capita income
                chg_lambda(r) lagrance multiplier
                chg_fd(i,r) final demand
                chg_absorption(i,r) absorption
                chg_production(i,r) production
                chg_Index(i,r) Price index P
                chg_supplier(i,r) Supplier effect S
                chg_fprice(f,r) factor price
*                chg_wage(f_Lab,r) wages
* only relevant for the sectors in res_sect
		chg_resprice(i,r)  price of sector specific resource
             
* intermediate nests to simplify notation (not necessary)
		res_nonres(i,r)   CES -- incl. income shock
		factorcost(i,r)	  C-D
		res_nonres_forfmc  sameas above but without income shock
 ;

equations eq_budget, eq_tradeprod, eq_index, eq_wage, eq_res_nonres(i,r), eq_factorcost , eq_FMC_nonres, eq_FMC_res, eq_res_nonres_forfmc ;


* Combining budget and demand:
* CLM - quad  - beware of sign of b (estimated with minus sign)
$if "%spec%"=="NH_CLMquad"	eq_budget(r)..       chg_income(r) =e= sum(i,dem_sh(i,r) *  
$if "%spec%"=="NH_CLMquad" 	exp(sigma * log(chg_income(r)) + (1-sigma) * log(chg_index(i,r)) +  (1-sigma) * 
$if "%spec%"=="NH_CLMquad" 	( rho(i) * log(chg_lambda(r)) - b(i) * ( ( log(chg_lambda(r) * U(r)))*( log(chg_lambda(r) * U(r))) - logUsquare(r))))
$if "%spec%"=="NH_CLMquad" 	);

*
* CLM - shifter 
$if "%spec%"=="NH_CLMshift"	eq_budget(r)..       chg_income(r) =e= sum(i,dem_sh(i,r) *  
$if "%spec%"=="NH_CLMshift" 	exp(sigma * log(chg_income(r)) + (1-sigma) * log(chg_index(i,r)) +  (1-sigma) * 
$if "%spec%"=="NH_CLMshift" 	( rho(i) * log( ( chg_lambda(r) * U(r) + b(i)) / (U(r) + b(i))))
$if "%spec%"=="NH_CLMshift" 	));

* where log g = ( rho(i) * log( ( chg_lambda(r) * U(r) + b(i)) / (U(r) + b(i))))


* CLM standard
$if "%spec%"=="NH_CLM"	eq_budget(r)..  chg_income(r) =e= sum(i,dem_sh(i,r) *  chg_income(r)**sigma * chg_lambda(r)**(rho(i))  * chg_index(i,r)**(1-sigma)  );

* TO DO
* if homothetic
$if "%spec%"=="H" eq_budget(r)..       chg_income(r) =e= sum(i,dem_sh(i,r) * chg_lambda(r)**(-1) * chg_index(i,r)**(1-sigma(i))  );

*eq_budget(r)..		chg_income(r) =e= sum(i,dem_sh(i,r) * chg_lambda(r)**(-sigma(i)) * chg_index(i,r)**(1-sigma(i))  );


* Combining Trade, production and absorption:
eq_tradeprod(i,r)..  chg_production(i,r) =e=  sum(s,prod_sh_adj(i,r,s) * chg_supplier(i,r) * bitc_shock(r,s)**(-theta(i)) *
                         chg_index(i,s)**theta(i) * (gosh_f(i,s)*chg_fd(i,s) + sum(j$gosh_m(i,j,s),gosh_m(i,j,s)*chg_production(j,s))) );


* Combining index and supplier effect:
*eq_index(i,s)..      chg_index(i,s)**theta(i) * sum(r,trade_sh(i,r,s) * bitc_shock(r,s)**(-theta(i))  *
*                     prod(f$beta(f,i,r),(chg_fprice(f,r)/income_shock)**beta(f,i,r))**(-theta(i)) * prod(j,chg_index(j,r)**intersh(j,i,r))**(-theta(i)) ) =e= 1;


eq_index(i,s)..      chg_index(i,s)**theta(i) * sum(r,trade_sh(i,r,s) * bitc_shock(r,s)**(-theta(i))  *
			prod(j,chg_index(j,r)**intersh(j,i,r))**(-theta(i)) *
			(res_nonres(i,r))**(-theta(i) * beta_all(i,r) ) 
			) =e= 1;

* middle nest (CES) CHI (in text)  


eq_res_nonres(i,r)..  res_nonres(i,r) =e=  (beta_Res(i,r) * ( (chg_resprice(i,r)* (1/income_shock))*income_shock_toresourcefactor
										+(chg_resprice(i,r))*(1- income_shock_toresourcefactor)
					 )**(1-nu(i,r)) + (1-beta_Res(i,r)) * (  factorcost(i,r) /income_shock   )**(1-nu(i,r)) )**(1/(1-nu(i,r)));


* income_shock is applied above
eq_factorcost(i,r) ..  factorcost(i,r) =e= prod(f$f_noRes(f),(chg_fprice(f,r))**beta_fnores(f,i,r));		
  

* !! cannot use the same nest for fmc's because we don't want the shock here..
* this is and intermdiate value to simplify notation
eq_res_nonres_forfmc(i,r)	 ..  res_nonres_forfmc(i,r) =e=  (beta_Res(i,r) * (chg_resprice(i,r) )**(1-nu(i,r)) + (1-beta_Res(i,r)) * (  factorcost(i,r) )**(1-nu(i,r)) )**(1/(1-nu(i,r)));



eq_FMC_nonres(f,r)$f_noRes(f) ..	chg_fprice(f,r) =e= sum(i$(not res_sect(i)),  fact_sh(f,i,r) * chg_production(i,r)) 
								+ sum(i$(res_sect(i)),  fact_sh(f,i,r) * chg_production(i,r)* 
								(factorcost(i,r)**(1-nu(i,r)))  * (res_nonres_forfmc(i,r)**(nu(i,r)-1))) ; 




eq_FMC_res(i,r)$res_sect(i) ..		chg_resprice(i,r)**nu(i,r) =e= chg_production(i,r) * (	res_nonres_forfmc(i,r)**(nu(i,r)-1));


* FOUR STEPS PER LOOP:
* 1) SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:
* 2) SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:
* 3) SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:
* 4) MANUALLY ADJUSTING FACTOR PRICES



* MODELS:

* 1) MODEL TO SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:
*model counterfactMCP_PE_justprices / eq_index.chg_index /;
model counterfactMCP_PE_justprices / eq_index.chg_index, eq_factorcost.factorcost, eq_res_nonres.res_nonres /;


counterfactMCP_PE_justprices.iterlim = 1000000;

* 2) MODEL SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:
model counterfactMCP_PE_demand / eq_budget.chg_lambda /;
counterfactMCP_PE_demand.iterlim = 1000000;

* 3) MODEL FOR SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:
model counterfactMCP_PE_fixedprice / eq_tradeprod.chg_production/;
counterfactMCP_PE_fixedprice.iterlim = 1000000;

* 4) MODEL FOR SOLVING FOR factor prices, TAKING PRODUCTION AS GIVEN:
model counterfactMCP_PE_wage / eq_FMC_nonres.chg_fprice, eq_FMC_res.chg_resprice,  eq_factorcost.factorcost, eq_res_nonres_forfmc.res_nonres_forfmc /;
counterfactMCP_PE_wage.iterlim = 1000000;


*------------

* APPROXIMATING STARTING VALUE FOR RESOURCE PRICE (speeds up convergence)

parameters  chg_fd_approx, chg_absorption_approx, chg_fprice_justfd, chg_fprice_approx,  log_Y_qty_approx;

chg_fd_approx(i,r) = 1 + (incelast("DIRECT",i,r) - 1) * (income_shock - 1);

display chg_fd_approx;


log_Y_qty_approx("approx",i,r ) =  log(income_shock) * (incelast("total",i,r))  ;
log_Y_qty_approx("approx",pe,r ) =  log(income_shock) * (incelast("total",pe,r)) * (1 - 1 /(1+ supply_elast(pe,r))) ;

display log_Y_qty_approx;


equation eq_FMC_res_approx;
eq_FMC_res_approx(i,r)$res_sect(i) ..		chg_resprice(i,r)**nu(i,r) =e= chg_production(i,r) * ( (beta_Res(i,r) * (
						chg_resprice(i,r)
						)**(1-nu(i,r)) + (1-beta_Res(i,r)) * (  factorcost(i,r)   )**(1-nu(i,r)) )**(-1));



model eq_FMC_res_presolve / eq_FMC_res_approx.chg_resprice /;


chg_production.FX(i,r) = ((income_shock-1) * supply_elast(i,r) + 1) /income_shock;

factorcost.FX(i,r) = 1;

res_nonres.L(i,r) =1;
chg_resprice.L(i,r) =1;


solve eq_FMC_res_presolve using MCP;

chg_production.UP(i,r) = 1000;
factorcost.UP(i,r) = 1000;

chg_production.LO(i,r) = 1/1000;
factorcost.LO(i,r) = 1/1000;



*------------------------



parameter count, Sqr_stat, Abs_stat;
count = 0;
Abs_stat = 100;
Sqr_stat = 100;

parameters      Sqr_stat_rep(T)
                Abs_stat_rep(T)
                chg_income_rep(T,r) per capita income
                chg_lambda_rep(T,r) lagrance multiplier
                chg_fd_rep(T,i,r) final demand
                chg_absorption_rep(T,i,r) absorption
                chg_production_rep(T,i,r) production
                chg_production_qty_rep(T,i,*) production
                chg_Index_rep(T,i,r) Price index P
                chg_supplier_rep(T,i,r) Supplier effect S
                chg_fprice_rep(T,f,r) factor price
                chg_resprice_rep(T,i,r) factor price
                chg_income_ini(r) per capita income
                chg_fprice_ini(f,r) factor price
           
                tradecosts_rep(r,s) trade costs
                chg_resprice_import(i,r) Storing resource factor price
                chg_fprice_import(f,r) Storing factor price;


tradecosts_rep(r,s) = 0;


chg_fprice_import(f,r) = 1;

chg_resprice_import(i,r) = chg_resprice.L(i,r);


bitc_shock(r,s) = tcost_shock;
bitc_shock(r,r) = 1;

tradecosts_rep(r,s)= bitc_shock(r,s);

display bitc_shock;


*$ontext
* Resetting count and convergence stats:
count = 0;
Abs_stat = 100;
Sqr_stat = 100;


*Bounds:
chg_income.LO(r) =  0.001;
chg_lambda.LO(r) =  0.001;
chg_fd.LO(i,r) =  0.001;
chg_absorption.LO(i,r) =  0.001;
chg_production.LO(i,r) =  0.001;
chg_Index.LO(i,r) =  0.001;
chg_supplier.LO(i,r) =  0.001;
chg_fprice.LO(f,r) =  0.001;

chg_resprice.LO(i,r) =  0.001;

res_nonres.LO(i,r) =  0.001;
res_nonres_forfmc.LO(i,r) =  0.001;


factorcost.LO(i,r) =  0.001;


chg_income.UP(r) = 1000;
chg_lambda.UP(r) = 1000;
chg_fd.UP(i,r) = 1000;
chg_absorption.UP(i,r) = 1000;
chg_production.UP(i,r) = 1000;
chg_Index.UP(i,r) = 1000;
chg_supplier.UP(i,r) = 1000;
chg_fprice.UP(f,r) = 1000;

chg_resprice.UP(i,r) = 1000;

res_nonres.UP(i,r) =  1000;
res_nonres_forfmc.UP(i,r) =  1000;


factorcost.UP(i,r) =  1000;


*Initial value:
chg_income.L(r) = 1;
chg_lambda.L(r) = 1;
chg_fd.L(i,r) = 1;
chg_absorption.L(i,r) = 1;
chg_production.L(i,r) = 1;
chg_Index.L(i,r) = 1;
chg_supplier.L(i,r) = 1;
chg_fprice.L(f,r) = 1;

chg_resprice.L(i,r) = 1;

res_nonres.L(i,r) =  1;
res_nonres_forfmc.L(i,r) =  1;


factorcost.L(i,r) =  1;


* INITIAL ADJUSTMENT FOR INCOME COUNTERFACTUAL:
* 
chg_Index.L(i,r) = 1/income_shock;


chg_fprice.FX(f,r) = 1;


chg_resprice.FX(i,r) = chg_resprice.L(i,r);

chg_Income.FX(r) = 1;

solve counterfactMCP_PE_justprices using MCP;


chg_index.FX(i,r) = chg_index.L(i,r);


chg_supplier.FX(i,r) =	 prod(j,chg_index.L(j,r)**intersh(j,i,r))**(-theta(i)) *
			(res_nonres.L(i,r))**(-theta(i) * beta_all(i,r) ) ;



display chg_index.L;

solve counterfactMCP_PE_demand using MCP ;


$if "%spec%"=="NH_CLMquad" chg_fd.FX(i,r) = exp(sigma * log(chg_income.L(r)) + (1-sigma) * log(chg_index.L(i,r)) +  (1-sigma) * 
$if "%spec%"=="NH_CLMquad" 	( rho(i) * log(chg_lambda.L(r)) - b(i) * ( ( log(chg_lambda.L(r) * U(r)))*( log(chg_lambda.L(r) * U(r))) - logUsquare(r))));

$if "%spec%"=="NH_CLMshift" chg_fd.FX(i,r) = exp(sigma * log(chg_income.L(r)) + (1-sigma) * log(chg_index.L(i,r)) +  (1-sigma) * 
$if "%spec%"=="NH_CLMshift" 	( rho(i) * log( ( chg_lambda.L(r) * U(r) + b(i)) / (U(r) + b(i)))));

* CLM standard
$if "%spec%"=="NH_CLM" chg_fd.FX(i,r) = chg_income.L(r)**sigma * chg_lambda.L(r)**(rho(i)) * chg_index.L(i,r)**(1-sigma);


*  (not used)
$if "%spec%"=="H" chg_fd.FX(i,r) = chg_income.L(r)**sigma * chg_lambda.L(r)**(rho(i))* chg_index.L(i,r)**(1-sigma(i));

chg_lambda.FX(r) = chg_lambda.L(r);

display chg_lambda.L, chg_fd.L;

solve counterfactMCP_PE_fixedprice using MCP ;
display chg_production.L;

chg_absorption.L(i,s) = gosh_f(i,s) *chg_fd.L(i,s) + sum(j$gosh_m(i,j,s),gosh_m(i,j,s) *chg_production.L(j,s));




* T LOOPS:
loop(T,

*for counterfactual in tau:
If ((Abs_stat > 0.0000001),


* 1) SOLVE FOR PRICES, TAKING FACTOR PRICES AS GIVEN:

chg_supplier.LO(i,r) = 0.001;
chg_supplier.UP(i,r) = 1000;
chg_index.LO(i,r) = 0.001;
chg_index.UP(i,r) = 1000;

solve counterfactMCP_PE_justprices using MCP ;


chg_index.FX(i,r) = chg_index.L(i,r);


chg_supplier.FX(i,r) =	 prod(j,chg_index.L(j,r)**intersh(j,i,r))**(-theta(i)) *
			(res_nonres.L(i,r))**(-theta(i) * beta_all(i,r) ) ;


* 2) SOLVING FOR DEMAND, TAKING PRICES AND INCOME AS GIVEN:

chg_fd.LO(i,r) = 0.001;
chg_fd.UP(i,r) = 1000;
chg_lambda.LO(r) = 0.001;
chg_lambda.UP(r) = 1000;

solve counterfactMCP_PE_demand using MCP ;


* clm flex
$if "%spec%"=="NH_CLMquad"	  chg_fd.FX(i,r) = exp(sigma * log(chg_income.L(r)) + (1-sigma) * log(chg_index.L(i,r)) +  (1-sigma) * 			    
$if "%spec%"=="NH_CLMquad"      ( rho(i) * log(chg_lambda.L(r)) - b(i) * ( ( log(chg_lambda.L(r) * U(r)))*( log(chg_lambda.L(r) * U(r))) - logUsquare(r))));

$if "%spec%"=="NH_CLMshift"   chg_fd.FX(i,r) = exp(sigma * log(chg_income.L(r)) + (1-sigma) * log(chg_index.L(i,r)) +  (1-sigma) * 
$if "%spec%"=="NH_CLMshift" 	( rho(i) * log( ( chg_lambda.L(r) * U(r) + b(i)) / (U(r) + b(i)))));

$if "%spec%"=="NH_CLM"       chg_fd.FX(i,r) = chg_income.L(r)**sigma * chg_lambda.L(r)**(rho(i)) * chg_index.L(i,r)**(1-sigma);

$if "%spec%"=="H" chg_fd.FX(i,r) = chg_income.L(r)**sigma * chg_lambda.L(r)**(rho(i))* chg_index.L(i,r)**(1-sigma(i));

chg_lambda.FX(r) = chg_lambda.L(r);


* 3) SOLVING FOR PRODUCTION, TAKING PRICES AS GIVEN:

solve counterfactMCP_PE_fixedprice using MCP ;

chg_absorption.L(i,s) = gosh_f(i,s) *chg_fd.L(i,s) + sum(j$gosh_m(i,j,s),gosh_m(i,j,s) *chg_production.L(j,s));


* 4) SOLVING FOR WAGES, TAKING PRODUCTION AS GIVEN:

chg_production.FX(i,r) = chg_production.L(i,r);

chg_fprice.LO(f,r)$f_noRes(f) = 0.001;
chg_resprice.LO(i,r)$res_sect(i) = 0.001;

chg_fprice.UP(f,r)$f_noRes(f) = 1000;
chg_resprice.UP(i,r)$res_sect(i) = 1000;


solve counterfactMCP_PE_wage using MCP ;

chg_production.LO(i,r) = 0.001;
chg_production.UP(i,r) = 1000;



* 4) CONVERGENCE STATS + MANUALLY ADJUSTING FACTOR PRICES:



Sqr_stat = ( sum((r,f)$f_noRes(f), sqr(1-  chg_fprice.L(f,r) /  chg_fprice_import(f,r)   )) /  sum((r,f)$f_noRes(f),1)
		+ sum((r,i)$((not no_ene_prod(i,r)) AND res_sect(i)), sqr(1-  chg_resprice.L(i,r) /  chg_resprice_import(i,r) ) ) /  sum((r,i)$((not no_ene_prod(i,r)) AND res_sect(i)),1) )**(1/2) ;


Abs_stat = sum((r,f)$f_noRes(f), abs(1-  chg_fprice.L(f,r) /  chg_fprice_import(f,r) ) ) /  sum((r,f)$f_noRes(f),1) 
		+ sum((r,i)$((not no_ene_prod(i,r)) AND res_sect(i)), abs(1-  chg_resprice.L(i,r) /  chg_resprice_import(i,r) ) ) /  sum((r,i)$((not no_ene_prod(i,r)) AND res_sect(i)),1)  ;


count = count + 1;

display Sqr_stat, Abs_stat, count, chg_fprice.L;




* weights can be 1/3 and 2/3rd for income counterfactual.
*chg_fprice.FX(f,r) = (1/3) * chg_fprice.L(f,r) + (2/3) * sum(i,fact_sh(f,i,r) * chg_production.L(i,r));


chg_fprice.FX(f,r) = (1/5) * chg_fprice.L(f,r) + (4/5) * chg_fprice_import(f,r);
*HERE FIRST TErM IS THE NEW ONE 


* 1/10 works well  (2/10 doesnt work)
chg_resprice.FX(i,r) = (1.5/10) * chg_resprice.L(i,r) + (8.5/10) * chg_resprice_import(i,r);

* compute change in income here

chg_income.FX(r) = sum(f$f_noRes(f),endow_sh(f,r)*chg_fprice.L(f,r))
	 + sum(i$res_sect(i),  endow_sh_res(i,r) * chg_resprice.L(i,r));


chg_fprice_import(f,r) = chg_fprice.l(f,r);
chg_resprice_import(i,r) = chg_resprice.l(i,r);

	

* Reporting:
Sqr_stat_rep(T) = Sqr_stat;
Abs_stat_rep(T) =  Abs_stat;
chg_lambda_rep(T,r) = chg_lambda.L(r);
chg_fd_rep(T,i,r) = chg_fd.L(i,r);
chg_absorption_rep(T,i,r) = chg_absorption.L(i,r);
chg_production_rep(T,i,r) = chg_production.L(i,r);
chg_Index_rep(T,i,r) = chg_Index.L(i,r);
chg_supplier_rep(T,i,r) = chg_supplier.L(i,r);

chg_fprice_rep(T,f,r)$f_noRes(f) = chg_fprice.L(f,r);

chg_resprice_rep(T,i,r)$res_sect(i) = chg_resprice.L(i,r);

chg_income_rep(T,r) = chg_income.L(r) ;

chg_production_qty_rep(T,i,r) = chg_production_rep(T,i,r) / chg_Index_rep(T,i,r);

chg_production_qty_rep(T,i,"world") = sum(r, production(i,r) * chg_production_rep(T,i,r) / chg_Index_rep(T,i,r)) / sum(r, production(i,r));


);
* end of "if" statement


display count;

);
* end of loop across T's


* dump all parameters 
execute_unload 'results%SLASH%CFresults_co2_%spec%_%demandest%_%prodspec%_fitteddem.gdx'





